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TECHNICAL NOTE 3511 


EXTRAPOLATION TECHNIQUES APPLIED TO MATRIX METHODS 
IN NEUTRON DIFFUSION PROBLEJffi 
By Robert R. McCready 


SUMMARY 

A general matrix method is developed for the solution of 
characteristic -value problems of the type arising in many physical appli- 
cations . The scheme employed is essentially that of Gauss and Seidel 
with appropriate modifications to mate it applicable to character istic- 
value problems . An iterative procedure produces a sequence of estimates 
to the answerj and extrapolation techniques, based upon previous be- 
havior of iterants, are utilized in speeding convergence. Theoretically 
sound limits are placed on the magnitude of the extrapolation that may be 
tolerated . 

This matrix method is applied to the problem of finding criticality 
and neutron fluxes in a nuclear reactor with control rods. The two- 
dimensional finite-difference approximation to the two-group neutron- 
diffusion equations is treated. Results for this example are indicated. 

The calculations 'vrere performed on the IBM card-programmed calciilator. 


INTROmCTIQN 

A general matrix method is developed for the solution of 
characteristic-value problems of a type arising in many physical appli- 
cations. The method of this paper is essentially that of Gauss and 
Seidel (ref. l), which itself is but a special case of the method of con- 
jugate gradients (ref. 2), The adaptation of the Gauss-Seidel technique 
to the characteristic -value problon calls for a means of coitputing suc- 
cessive estimates of the characteristic value eis well as the vector. 

This calculation is made to rely upon Turner’s technique (ref. 3) for 
assigning a meaning to the ratio of two vectors. 

Extrapolation techniques are also employed to speed up the conver- 
gence of the Iterative process. One of these is based on Turner's orig- 
inal formula (ref. 3), and the other is a sli^tly more conplicated 
modification. 
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The number of iterations reqtilred for convergence is not studied 
theoretically here as in the "n-step" methods^ hut the minimization of 
a suitable form at each step is derived. 

The method is applied to two-group neutron-diffusion equations. 

The calculations were performed on IBM equipment at the MCA Lewis 
laboratory by Mary J. Winter. 


SYMBOLS 

The following symbols are used in this report: 
A^B^L^U matrices 

axial leafcage 
D,E,F,G^J,X vectors 


^th 

2 


h 

LtJ 

N 

^th 


sgn u = 


w_,co 


grid dimension 
indices 

thermal multiplication constant 

average square slowing down length for fast neutrons 

average square diffusion length for thermal neutrons 

number of nuclei per cc 

resonance escape probability 

radial coordinate 

core radius 

' 1 u > 0 
-1 u < 0 
u = 0 

reflector thickness 
wei^t functions 
characteristic value 
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k+1 


X 


T 


deviation at i^^^ point of kth iteration (eq. (80)) 
difference 

actual dan^jing rate 
bulk dartq)lng rate 
neutron fluxes 


Parameter groupings: 


Subscripts ; 
f faat 


% 

, Hr,tbo 1 


1 -n 2 

0 = + 

Wo_ 

\r,thg ®**>0 2 


I 2 ®z 


Lf, 


g = 


2-^®z 


" " Lf = 






tb tbermal 
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tr transport 

0 reactor 

1 reflector 

2 rod 


THE METHOD 

Matrix fo rmu lation . - Consider tlie matrix equation 

AX = rBX (l) 

■where A and B are n x n matrices, X is an n-component vector, and 
the characteristic value y is a scalar to he determined. A may he 
separated into the sum of two triangular matrices L and U where L 
contains all the diagonal elements of the original matrix A. 

This separation, which anticipates the Gauss-Seidel process, is 
effected in the following manner: 



A = L + U 


(2) 

hj = "ij 


J >i 

( 3 ) 

"ij '“u 

J >ij ° 

j i 

( 4 ) 


If L is a nonsingular matrix (always true if ^ 0 for al 1 

l), equation (l), modified to 

(L + U)X = rBX (5) 

may he multiplied hy L“^, giving 

(I + L"^)X = rL"^BX (6) 

For a given X, the quantities L~^UX and L“^BX of equation (6) 
may he ceilculated without the actual formation of L“^. This fact, idiich 
is very helpful for systems containing large matrices, arises in the fol- 
loi-fing manner and depends upon the triangular nature of L. Let D he 
the vector defined hy 


D = L"^ 


( 7 ) 
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Then 


LD = UX =■ C 


( 8 ) 


whence 


^11*^1 “ ‘^l 


( 9 ) 


gives since all the Cj^ can he cong>uted from U and X in equa- 

tion (8). Then 


^2l'^l ^22*^ “ (lO) 

gives d 2 , and 

^3A ^32*^ ^33*^ “ *^3 

gives dj, and so forth, so that L“^ need not be con 5 >uted in order to 
obtain L"^UX. The same argument applies to L“%X. 

Iterative scheme . - Eqtiation (6) may be v/ritten 


X = rL"^X - L’^ 

which may be interpreted aa defining the iterative scheme 

Vi “ - l'H 


( 12 ) 


(13) 


in which is an estimate to y that can be calculated from Xj^. 

To obtain form the inner product of the vector sgn L~^BX^ with 

each side of eqimtion (6)j thus. 


k+1 


(sgn (I + L~%)X)^) 

(sgn L-^BXj^, L"^X^) 


(14) 


Equations (l3) and (l4) axe the basic equations of the iterative scheme. 
Given any Xj^, is computed from equation (l4) and aai 

are planed in (l3) to yield the next Iterant This process is 

repeated until X^, and converge. 


Some normalization is necessary in problems of a homogeneous nature. 
The sln 5 )lest method of normalization is to adjust a permanently specified 
coordinate of X^^. to unity before beginning each iteration. This is 
acccm 5 )lished by dividing each element of the vector by the specified 
coordinate . 
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IHie ratio defined by equation (l4) was chosen for sin 5 )licity of 
calculation on available punched-card equipment. That ratio can be com- 
pared to the Raylei^ quotient (for eq. (l3)) 


'k+1 J^) 


(15) 


where 


Jk = i-'X 


(16) 


Gjj. = (I + L“%)Xjj. 


(17) 


by noting that each of the relations (l4) and (l5) constitutes a weighted 

, N (i)' 

s-um of local ^point by point; values of ^k+1* These local values 

are defined by 



_(1) = Ji)/,(i) 


= r 


k+l 


= 


(18) 


where g^^^ and are the ith components of and 

tively. The weighted average associated T^rith (l5) is 


Jjj., respec 



where 



while that associated with (l4) is 


I'k+l - E “i ^k+i 

i 


(19) 


(20) 


( 21 ) 
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where 


“i 




E 

m 






( 22 ) 


Equation (l5) selects that value of which m i nimizes the sum 

of the squares of the residuals of equation (6) when that quantity is 
thought of as a function of Til simi of the squares of the res id- 

iCrX 

uals is not, of course, the only quadratic form that is suitable for 
minimization (ref. 2). Consider the expression 






,W 


LL‘Jk 


sgu j 




1/2 




1/2 


(23) 


i/hich is zero for Xj^. = X and = y. This generally positive quan- 

tity can he minimized (made closer to its ultimate value zero) by setting 


r 


(sgn Jj^, Gj^) 


k+1 (sgn J, 


k^ 




(24) 


which is equation (l4) expressed in terms of J and G. 

Example . - To illiastrate the convergence of this method in a spe- 
ci£Ll case, consider the problem of equation (l) with 



(25) 


and 


B = 


which has the real solution 

v(3) 


10 0 0 
^^0 0 O 0 ^ 

X^^^ = 1.020070, X^^^ = 1.329658, 


(26) 


1.000000, x(^^ = 1.109886; r = 0.549429 and two solutions >rith 
conplex characteristic values . This solution was f oixnd by the ordinary 
process of solving the characteristic equation. 
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This problem, 'vjas solved ±a 15 iterations starting with an initial 
guess of Xq = (l0_, 100,1, 1000 ) . The values of successive iterants, 
together with those of y, are listed in the following table. The iter- 
ants are normalized so that =1 at the start of each iteration: 

k 


k 

41) 

42) 

44) 

"'^k+l 

0 

10 

100 

1000 


1 

-0.692996 

1.069635 

1.023212 

-10.634589 

2 

.271412 

.438315 

.812771 

-1.352353 

3 

.848870 

1.565152 

1.188384 

-.097105 

4 

1.132470 

1.360721 

1.120240 

.526554 

5 

1.006120 

1.292707 

1.097568 

.592807 

6 

1.010639 

1.337044 

1.112348 

.538366 

7 

1.024544 

1.331918 

1.110639 

.547271 

8 

1.019872 

1.328104 

1.109368 

.551383 

9 

1.019593 

1.329848 

1.109949 

.549107 

10 

1.020236 

1.329788 

1.109929 

.549291 

11 

1.020078 

1.329596 

1.109865 

.549510 

12 

1.020048 

1.329660 

1.109886 

.549423 

13 

1.020076 

1.329665 

1.109889 

.549422 

14 

1.020071 

1.329656 

1.109885 

.549432 

15 

1.020070 

1.329658 

1.109886 

.549430 


Extrapolation technique . - If, instead of four components, the iter- 
ant vector has maaay ccm^jonents, techniques of extrapolation are usually 
desirable to speed convergence of the process. The technique enqoloyed 
here, which is due to Turner (ref. 3), attempts to evaluate a bulk damp- 
ing rate 'vjhich describes, in an average way, the over-all trend of the 
individvial components of the iterant vectors. 


Assume that each iterant Xj^. is made up of the sum of the solution 
X and two error vectors and Fjj. satisfying the damping relations 

C27) 


and 


•k+1 


= - tE 


(28) 
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Then the following relations hold: 


w 

05 

05 


W 

H 

O 


One may ccng)ute 


^0 = 

X + Eq + Fq 

(29) 

Xi = 

X+'CEq -tFq 

(30) 

Xg = 

X + + t%-q 

(31) 

^3 = 

X + t\ - 

(32) 

t2 = 

'\-^o 

(33) 


The "vector division" implied in f33) is possible because^ under the 
initial assim^jtion of error behavior ((27) and (28)), the vectors 
X5 - Xg and X-j^ - Xq are coUineax and therefore differ only in length. 


If the error vectors are eliminated from equations (30) and (32), 
one obtain^ 


X = 



(34) 


■vdiich gives the answer as a linear combination of the alternate iterants 

X.| a.Tifi ^3 * 


The preceding analysis sijggests that a formula analogous to (34) be 
lised to estimate the answer. The difficulty here is that equation (33) 
may be meaningless when equations (29) ito (32) do not hold. To circumvent 
this difficulty, a method of con5)utlng Is needed. Tov/ard this end, 

define 8^'^} by means of 


(i) ^ j 5-(<3) vCj) 


^+i = 


k+1 


■O' » ^ 
■ ^ 


(35) 


and define 


by means of 


X 


2 



(36) 
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0?he direct analogy to equation (l4) will "be noticed. Eqixation (36) per- 
mits con^nitation, in an average t'ray, of the danrpiag of the error vectors. 
With available^ the extrapolated value X' of X is computed from 


X’ 



(37) 


In case the error is dan 5 )ing exactly as assumed in (27) and (28), equa- 
tion (36) gives the value indicated hy (33) and equation (37) reduces to 
(34:)j that is, X' becomes the answer X. 


Since the ideal dacping behavior is rarely an actuality, it is of 
interest to examine the effect of the preceding process on error con^on- 
ents. Suppose that Xq is more adequately represented by 

n 

^ ^ 

i=l 


\7here 



has a (damp ing rate (positive or negative) of 
X 3 = X + 

i=l 



Then 


(39) 


and 


Xj_ = X + 


E 


X 


(4=0) 


hold. The extrapolation indicated in equation (37) now yields the fol- 
lowing relation between the estimate X' and the answer X: 


X' 


= X + 


i=l 1 - ^ 


(41) 


This interpretation is useful, since it indicates the damping effect 
ijpon the errors of three iterations and one extrapolation. 




If, for simplicity, one of the errors E^^^ and its damping rate 
are designated by E and X, respectively, then 


R(X,t) 


(X^ - T^)X^'^ 
1 -t2 


(42) 
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gives the dau^jing of this error component as a resiilt of j iterations 
and one extrapolation. The "eidjreme" value of E (actually that value 
farthest from zero; i.e.^ farthest from maximum damping) may be found by 
setting 


dR 
dX = 


.jX-^ ^ 

1 -t2 


= 0 


(43) 


This yields 



(44) 


as the equation to be solved for the values of X which are associated 
•vdLth the errors that receive the minimum damping frcm the process of j 
iterations and one extrapolation. Equations (42) and (44) give 
the extreme value of R, as a function only of r and j: 


I^extC'^) 


1 fi - 2y /^ 

l_T^‘5"2y j J 


(45) 


To find the value of so that this function cannot exceed the 

boimds dbl - that is, so that the slowest dancing ccsi 5 )onent (and hence 
all components) cannot increase through extrapolation - T must be less 
in absolute value than the least of the roots of 



(46) 


p 

If only such T are used, the convergence of the process cannot be 
in^jaired by the extrapolation. 


Suppose now that the previous value of 
formula 


R(X,t) is replaced by the 


r(x^) 



(47) 


In formula (47), j ^ 4. The second factor places a zero (max. 
damping) at just the points of mlnlnn^m dang)ing, that is, at the values 
of X determined by (44). If now dR/dX is taken as zero and the limit 

±1 is placed upon the resulting Rgxt(''^^^ limiting safe values of 

p 

r are obtained by finding the least of the roots of 
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|jr‘’ - 2(j - l)r'5’2 + U - 2)*5"^Jt«J ± [j - 2(j - l)^ + (j - 4)r^ = 0 

(48) 


■where r satisfies 


.2 ^ (3 - l)(j - 2 ) ± Vsjg - 12j + 4 

.2 


(49) 


The re‘\rised formula (47) has "both the effect of ensuring that no compo- 
nent -will he impaired in its daTi 5 )ing hy the extrapolation and also that 
the least rapidly darrping ccmponent recei-ves a zero contribution in the 
extrapolation . 


Since, as before, for some error ccmrponent E, 

Xj <= X + X‘5 e (so) 

Xj_2 = X + (51) 

X, , = X + X*^-^ (52) 

in -which X represents the answer, the specification of (47) as a damp- 
ing formula implies 


Y. - X + - 2(3 - + (J - 2)t^X*^'^ 

3 - 2(3 - 4 (3 - 


(53) 


where X’ is -bhe extrapolated value of Xj. If X*^E, and X'^"Sl 

are eliminated from (53) using the relations (50), (Sl), and (52), then 


3X. - 2(3 - 1)t2x. 2 + (j - 2 )t^._ 4 , ^ 

X' = ^ (54) 

3 - 2(3 - 1)t2 + (3 - 2)t^ 

Coanparison of (42) (-with 3 = 5} -with (37) on one hand and of (47) 
•with ( 54 ) on the other hand leads to the following valid rule of thumb 
to obtain the extrapolated -vulue of X for a given (damping f;mction: 

Replace the power X^ of X in -the daniping function by X^j the result- 
ing linear combination of alternate iterants is the formula for the ex- 
trapolated X. It is easily verified that the validity of this arises 
from the manner in which the error vectors are assumed to behave. 
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The smallest roots of equations (46) and (48) are listed in the 
follovring tables: 


Eq 

. (48) 

J 


4 

0.6667 

6 

.8745 

8 

.9213 

10 

.9427 


Eq 

. (46) 

J 


4 

0.8284 

6 

.8941 

8 

.9233 

10 

.9399 


These are the ijpper liirn'ts of the "safe" values of within the frame- 

work of the definition. 


APPLICATION TO REACTOR THEORY 

General remarks . - Multigroup reactor equations can he solved, in 
principle, hy the present method. The nimiber of coairponents in the vector 
solution, to be disciassed in detail later, is approximately equal to the 
product of the number of grid points and the nvimber of groups in the ntul- 
tigroup scheme. An extreme increase in the number of these elements 
lengthens the problem considerably. The calculations here are performed 
in accordance with two-group neutron-diffusion theory. 

The two-dimensional reactor with control rods, which is considered 
later, is siilted to two-group calculatibns, since the control rods are 
particularly effective on the thermal group and t>ro-group calcxilations 
are good for thermal assemblies . 

The foUoT'/lng illustration is introduced to show the general prin- 
ciples of the matrix setup in detail. These principles do not change 
for the more cong)licated two-dimensional problem that is treated later. 

A relatively simple one -dimensional problem has been chosen to illustrate 
the detailed setup and the consequent matrix. 

Example of tvro- group diffusion equations . - The one-dimensional dif- 
fusion equations for a reflected theimial reactor of slab geometry are 
(ref. 4) 

d^<p^ 

-a'<P„ = - rb'P^, 

X ’til 



( 55 ) 
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and 


for the core, and 


and 




th 


dx^ 


^"^th + d'Pf = 0 


d^q) 

-TT - 0 

dx 


d^. 


th 


dx 


S‘\h + = 0 


(56) 


(57) 


(58) 


for the reflector. 

The parameters a', h, c*, d, f, g', and m are defined in the list 
of symbols} r is the characteristic value of the system and equals 1 for 
criticality. Men r converges to a value other than unity, the uranium 
concentration is adjusted and the process repeated. 


The differential equations (55) to (58) are replaced hy finite- 

difference equations} the operation d <p/dx is estimated hy means of 
the approximate formula 


2 d^q>l 

1 7 - 

dx'^ 


0+1 0-1 j 


|Xi=X, 


(59) 


where the points of the region are numbered in order as grid points of 
a finite-difference net, and h is the distance between successive 
points. In the following, r^, is the core radiiis, the coniplete 

reactor radius, and point 6 lies on the interface: 


0 



0 12 
0 


3 4 5 6 7 


r 


c 


8 9 


r 


c+t 


The boundary conditions are that the fast and thermal flux have zero 
current across the plane of symmetry (x=0) . This condition 
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can be approximated by 



(60a) 


for both and The condition of continuity of currents at the 

Interface is met by approximating the derivatives in the expression 


, d«P / 
- \ro dj 



dx 


(r, + 0) 


(61) 


for both the fast and thermal fluxes. The remaining condition is that 
the flux be zero at the outer boundary. If the fast flux is designated 
by <P and the theimal by ijr, the system, becomes: 


Eq-uation 

J 

0 


1-5 


6 

7 


8 


<Pl - 




0 


= 0 


‘Plu.T + ‘P-1 1 - 2<P, 

J+1 J-1 J , ^ , 

2 ^ ‘PJ = - 


^ “Pg ~ ‘P 5 

Hr^fO 5 “ Hr,fl h 


‘P 7 -*P 6 


?8 + <*>6 - 2 ‘ P ? 


f’q>Y = 0 


‘P7 “ 2‘Po 

— "^'‘Pg “ ^ whicb incorporates cpg = 0 

h^ 


J 


(6?,) 


for the fast-balance equations and, for the thermals: 
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Equation 

i 

0 


1-5 


- +0 = 0 

'l' 1 -n + ♦i -1 ■ 1 

_i±i |-i cVj^ + d(Pj^ = 0 


h 

^6 - % 


% - '^p 


^tr,thO tL Hr^thl h 




> (63) 


■^8 + ^6 " ^ilry 


- g'iI/7 + nHp7 = 0 


- g'i|f + nKPo = 0 which incorporates 
80 




The variahles <Pq to <Pq and to may now he written as 

to X 3 and Xg to X^^, respectively. The matrix formulation of these 
equations is presented in figure 1. The following symbols have been 
introduced; 


It = 

(64) 

Aj^ => 2 /h^ + a 

(65) 

= 2 /h^ + f 

( 66 ) 

= 2 /h^ + c 

(67) 

Sjj = z/i? + g 

( 68 ) 

Hr,fl 

(69) 

„ Hr,thO 

” \r,thl 

(70) 


Recanitxilation. - To review the general application of the method^ 
to t^ ro-groT-ip reac^ equations, consider the following broad outline of 

this process: 

( 1 ) Write the two-group equations with the parameter T Introduced 
as a multiplier of the production term of the fast-balance equations. 
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(2) Express the differential equations hy their finite-difference 
approximations so that they become a linear algehrafc set of the type 
associated with equation (l) . 

(3) Perform, such iterations and extrapolations as necessary to ob- 
tain well-converged values of y an<i X. 

(4) Adjust the uranlTmi concentration and repeat step (3) using the 
original answer from (3) for the initial guess Xq. The concentration 
should be changed so that y 1. 

(5) Repeat (3) and (4) until y converges. If criticality is 
desired, change the concentration so that the converged values of y-»-l. 


TWO-DIMENSIORAL REACTOR WITH COBTROL ROIB 

Geometry of reactor . - The reactor (see fig. 2) is cylindrical and 
water-reflected with core congosed of aluminum, water, and uranium, which 
are assumed to be homogeneously mixed. The height of the reactor is 70 
centimeters, the outside radius 50 centimeters, and the core radius 32 
centimeters . 'Nine cadmium control rods are inserted in the core; one, a 
cylindrical rod of 2-cent±meter radius, is centered along the axis of the 
reactor. The remaining el^t rods are equally spaced on a radius of 24 
centimeters and are shaped so as to be bounded by coordinate stirfaces. 

Each of these rods extends over a radial distance of 4 centimeters and 
siibtends a central angle of 9°. 

The symmetry of this assembly is an inqportant factor in making solu- 
tion of the reeictor problem feasible. The flux in the 45° sector indi- 
cated in figure 2 is adequate to represent the flux in the entire reactor 
in fact, additional symmetry within the sector implies that only half 
the sector need be considered. The three-dimensional problem is made 
two-dimensional (confutation-wise) by estimating the neutron leakage in 
the axial direction due to the finite hei^t of the reactor. This is 
based upon an axial cosine distribution similar to the bare pile solution 
(eq. (75)). 

Composition and nuclear parameters . - The core volume is proportioned 
between the water (p = 1 g/cc) and aluminum by assuming a volume ratio of 
aluminum to water of 0.76. The nuclear diffusion constants for the core 
and reflector are listed in the following tables. The subscripts 0, 1, 2 
refer to the core, reflector, and rod regions of the reactor, respectively: 


• 
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Thermal 

Zone 

L 2 
^th 

\h 

^th 

\h 

0 

3.471 

0.815 

0.95 

1.675 

1 

8.3 

.43 

.98 



Fast 

Zone 


^f 

0 

64 

3.78 

1 

33 

3.43 

2 

-- 

4.35 


Parameters for the rod regions are unnecessary because of the sim- 
plified treatment of the rod, in ^diich the thermal neutron flux is assmed 
to vanish on the rod boundary and the radial and axial leakages are as- 
sumed to balance in the absence of fast neutron absolution processes . 

The thermal parameters in the preceding table are those associated 'vrf.th 

an atom ratio of — of 350; these, of course, change for different 

uranium concentrations. 


Equations and boimdary conditions . - The two-group equations (ref. 
4) for the core are taken to be 


7^ 


fO 


^ 1 ^ T. 2 V _ ^ Hr,thO J 

,2 ®2j<PfO " ^ X ^thO 2 -"thO 


= 0 (71) 


and 


V <P. 


thO 


+ (p._ + 


^tr.fO 


z/ thO X 


"thO 


tr,thO 


^thO ,2 ''’fO 
^fO 


(72) 


All the parameters of equations ( 71 ) and (72) are the ordinary nuclear 
ones, except the arbitrarily inserted y, which is a measure of the crit- 
icality and is equal to 1 for a critical assembly. 


In the reflector the two-group diffusion eqviations tale the foim 


7 tp 


fl 


■"fl ^ 


( 73 ) 


7^<P 


thl 


+ B \<p, , T + , 
z r^thl X. 


thl 


tr,fl 
tr , thl 


^thl ^2 ‘Pfl 


= 0 


fl 


( 74 ) 
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^/hile the fast -diffiis ion equation for the rods is taken to be 

Any change in rod boundary conditions would not affect the general princi- 
ples of the numerical scheme. As will be seen from the boundary condi- 
tions^ the thermal neutrons do not require a diffusion equation '>-rithin 

the rods. The region considered in the problem is a 22— sector (of the 

2 

circle or fig. 2)j one side of which, passes through the center of one of 
the outlying rods . This is illustrated in figure 3 . The symmetry of the 
over-all reactor in: 5 )lies that the normal derivatives of the flinc across 
the surfaces A and B are zero. This implies that the flux at all 
points of the circle of figure 2 can be found by solution for the flux 
only in the sector indicated in figure 3. The condition of -continuity 
of flinces and currents is Involved at the core-reflector interface^ in- 
dicated by C in figure 3. The vanishing of the fast and thermal flux 
at the outer boundary (d) is also required. The thermal flux is tal^en 
as zero on the edge of the control rod^ and the continuity of the fast 
flue and current is considered to hold on the core-rod interfaces . The 
details of the mathematical formulation of these conditions are deferred 
uatil the general discussion of the difference equations. 

Finite -difference equations . - In order to write the reactor equa- 
tions as finite -difference approximations,, the sector of figure 3 is 
divided iato a grid net of points . The flux is determined by solution 
of the linear algebraic system of equations I'dilch result from \vriting 
the finite -difference approximation to the fast- and thermal -diffusion 
eqxoations at each point. The grid arrangement iised in the present prob- 
lems is indicated In figure 4. The thermal flux (components 1 to 139 of 
the vector solution) has the foUo^rtng breakdovm into groups of compo- 
nents: reflector (l to 73), core-reflector interface (74 to 79), core 

(80 to 139). The fast flux (l40 to 29l) has the folloi'dng breakdown: 
reflector (l40 to 212), core -reflector interface (213 to 218), core (219 
to 284), control rods (231, 232, 237, 238, 243, 244, and 285 to 29l) . 

The number of components associated irtth the thermal and fast f luces, 
respectively, differs because of the assignment of boundary conditions at 
the control rods, vdiich brings the fast flux into a larger area of 
definition. 


V^(p 


For two-dimensional cylindrical geemetry, the operation of the form 


is given by 




( 76 ) 
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This form, is to be replaced by a difference operation that relates each 
point to its four nearest neighbors. If the point in question is desig- 
nated by the subscript zero and the others are 

1 


4 ' hg 0* hg • 2 r I 6 — 

3 

where h^ and hg are the grid widths in the r and 9 directions, 
respectively, then at ^ = <Pq the following approximation is used: 



With this designation (and barring certain exceptional points to be 
disctxssed later), one may move from point to point on the grid and wite 
equations of neutron balance for each of the two neutron groi; 5 )s. 

The following equations may be taken as typical illustrations ; 

Thermal -balance equation 93 (see fig. 4); 






^2 2rh^j 97 ^ ^Z^^2 '•'^92 -^94^ 


2 2 
+ + 


^hr^ hg2p2 


PthO ^234 = 0 
hO / tr,thO ^0 
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Fast-balance equation 234: 



In contradistinction to equations (78) and (79), there are certain 
special equations ■vdiich hold at the exceptional points referred to 
earlier. These equations result from one or more of the following 
conditions: 

(1) Continuity of currents at interfaces 

(2) Zero flux at the outside boundary 

(3) Zero current across planes of symmetry 

(4) Change in grid dimensions 

Condition (l) is treated by matching a suitable ratio of normal 
derivatives from either side of the interface. Each of these de- 
rivatives is evaluated by a five-point differentiation formula. Con- 
dition (2) is treated by witing the difference approximation to the 
diffusion equation for points adjacent to the outside boxmdary with 
zero replacing the fl\nc at the botmdeiry point. 

Condition (3) is accounted for by Vriting the diffusion eqiiation of 
a point on the plane assuming the same flux at grid points on either side 
of the plane. For condition (4), 0-wise interpolation foimulas are used 
to define fluxes at the points marked "X" on figure 4, and these are uti- 
lized, where needed, in siting diffusion equations in the finer net. If 
each equation of the set is written in order and the production terms are 
isolated as Illustrated in eqiiation (7l), then the rntrlx equation con- 
structed from the approximate finite difference may be written in the 
form of equation (l) . 

The matrix B is singular, largely consisting of zero elements with 
an essentially diagonal group of nonzero terms somewhat off the leading 
diagonal. The matrix A has a substantial number of nonzero elements 
crowded quite close to the leading diagonal. This latter situation is 
numerically desirable, as elements far from the leading diagonal tend to 
slow the convergence of numerical processes. 
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If criticality is desired, the concentration of fissionable material 
is adjtisted, after y and X have converged, and a tdiole new set of 
calculations is run mtil a new value for y is reached. This process 
may be continued until y =■ 1. 

The method can also be used to congoute reactivity changes; the cal- 
culation time is again shortened considerably if flux distributions are 
not demanded. 

Solutions of two-dimensional problem. - The results of the f.n.l mil h - 
tions of the supercritical (y = 0.948) case are shown in figures 5 to 7. 
Figure 5 gives the fact flux as a fronction of r for 0 = 0°, 9°, and 
18°. The control rods have no siibstantial effect on the fast fl^nc. Fig- 
\rre 6 gives the corresponding thermal flux and shows the localized effect 
of the control rods . Figure 7 presents iso-flux contours of the thermal 
fl'ux. The 0.19 contours in the reflector and the 0.234 contour in the 
core represent relative mayimumR . 


COMMEIKTS ON AEPIICATION OF THE METHOD 


A number of numerical quantities may be examined in an attempt to 
evaluate the degree of convergence of a system. One of the most natural 
of these quantities is the s\jm of the squares of the residuals. Another 
may be formed by considering the fact that, as the limit is approached, 

the ratio T must tend toward unity. This means that the devia- 

k+l/ k+1 

tion defined by 

.(i) 




k+l 


r. 


(80) 


h+l 


must tend tovrard zero. The average absolute value of the deviation, 
summed over a1 ^ points of the reactor is 



( 81 ) 


where n is the number of reactor points. 
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An illustration of the hehavior of this quantity is given in the 
foUovring table: 


k+1 


6 

15 

22 

31 

40 

49 




0.02199 

.02551 

.00254 

.00130 

.000793 

.00028 


k+1 


max 


0.2133 

.3156 

.0147 

.0113 

.0061 

.0027 


The iterants listed are those which just precede the extrapolation proc- 
ess. These are chosen so as to minimize the effect of fluctuations in- 
troduced by the extrapolation technique. 

These lllixstrative values come from the second general process; 
that iSj after y had converged to 1.2064, the concentration (and hence 
elements of the matrices A and B) was changed and a new series of itera- 
tions begun. This converged (more rapidly than the first run) to a 
vetlue of 0.948. 

To estimate the value of uranivmx concentration needed for the new 
run, the equation 


1.2064 kj-j^ (ol<i) = (new) 


(82) 


was used to compute a new k^.j^ from -vdiich to obtain a new concentration. 
This formula is em approximation, since the influence of a change in con- 

2 

centration upon is appreciable. The better rate of convergence of 

the second run is caused by the fact that the flux is relatively inde- 
pendent of the characteristic value so that the initial estimate for the 
second run was a relatively good one. 


The quantity 



than that of the vector 


reflects the convergence of 

X. 


y, which is faster 


In order to determine the degree of convergence of X, consider the 
quantities 




;+l 


■ E K;i| ■ E 

i=l i=l 


.(l) v(i) 


k+1 


(83) 
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j.(l) _ *^k+l 

k+1 n 


(84) 


and the mayimmi. 


8 . 


(i) 

'k+1 


designated hy 


quantities ane as follows: 



Typical values of these 


r = 1.206 

k+1 

\+l 


g(i) 

max 

g(i) 

k+1 

in 

122 

133 

144 

155 

0.0099 

.0083 

-.0049 

.0028 

.0022 

0.000034 

.000028 

.000017 

.000010 

.000008 

0.000269 

.000176 

.000157 

.0000904 

.0000709 


X = 0.948 

k+1 

Vi 


(i) 

k+1 


k+1 

max 

6 

15 

22 

31 

40 



0.3299 

.2139 

.0223 

.0044 

.0020 

o.oon37 

.000735 

.000077 

.000015 

.000007 

0.00918 

.00530 

.000204 

.0000666 

.0000245 


The sums of the squares of the residuals for the two cases y = 1.206 
and X - 0.948 axe as follows: 


r 

= 0.948 

k+1 


6 

i.53no“^ 

15 

3.51X10-4: 

22 

1.80X10"^ 

31 

1.22X10-7 

40 

1.99x10"^ 


r = 

1.206 

k+1 

2%^ 

in 

9.44x10"7 

122 

4.98x10"7 

133 

2.66X10-7 

144 

8.60X10"® 

155 

5.37X10"® 


Several general observations can be made about the process: 

(1) The n\nriber of iterations in this problem starting from an initial 
guess to a well -converged value of X was about 150 to 175. 

(2) In general, 8 to 10 iterations between extrapolations seem desir- 
able, as the use of too few iterations does not allow the establishment of 
a fairly uniform damping rate. 

(3) The extrapolation formula of equation (37) seems best for rou^ 
estimates ^diere error components are being damped rapidly; that of eqm- 
tion (54) seems to be superior for later extrapolations where one is 
closer to the solution. 


3Ld/ 
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(4) \Jhen ccai 5 )uted values of exceed the ijpper limits they may 

he replaced by the limit from the tables giving roots of equations (46) 
and (48) and extrapolation carried out^ or two more iterations perform^ 
with recon 5 )uted until it falls within prescribed limits. 

The following table gives the simi of the squares of the residuals 
for (a) direct iteration from "to (b) 8 iterations from 

oj followed by extrapolation with safe" -tdien "'r^ conputed" was too 

^ large, then iteration to X^g, (c) 8 iterations from ^31 foUoTOd by 

2 iterations and a test lantil ccmputed" was less than safe, " 

then extrapolation followed by iteration to X. 

4:y 


H 

o 


Case 


a 

l.SbxiO""^ 

b 

3 . 44 x 10 “^ 

c 

1.68X10“® 


Lewis Flight Prcpiilsion Laboratory 

National Advisory Committee for Aeronautics 
Cleveland, Ohio, May 12, 1955 
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Figure 1. - Matrix formulation of equations (62) and (63). 
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Figure 5. - Fact neutron flux for azimuth angle 6 = 0°, 9®, 18°. 
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